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SUMMARY 

An  algorithm  possessing  second-order  spatial  accuracy  has  been  presented  for 
solving  the  laminar  Incompressible  two-dimensional  Navier-St.okes  equations  In 
streamfunctlon-vortlclty  variables.  The  method  Is  either  second  order  accurate  in 
time  for  unsteady  flows  or  corresponds  to  Newton's  method  of  iteration  in  the  steady 
case.  The  formulation  Is  capable  of  handling  arbitrary  geometric  configurations  and 
the  alaorithm  contains  no  adjustable  parameters  nr  explictly  added  artificial 
dissipation.  All  terms  of  the  governing  equations  are  treated  implicitly  without 
incurring  any  factorization  error  as  is  commonly  generated  with  alternating- 
direction-implicit  schemes.  This  not  only  allows  Newton's  method  to  be 
employed  for  steady  calculations,  but  also  permits  implicit  differencing  of 
cross-derivative  terms  which  would  arise  for  nonorthogonal  grid  systems.  For  the 
unsteady  equations,  an  analysis  of  a  corresponding  constant  coefficient  system  is 
provided  which  indicates  numerical  stability.  Use  of  a  direct  solution  of  the 
linearized  equations  makes  the  method  noniterative  at  each  time  step  and  precludes 
the  appearance  of  spurious  time-like  terms  or  the  possibility  of  a  nonconverqing 
iteration  as  the  mesh  step  size  is  reduced.  In  addition,  machine  accuracy  is  not 
only  achieved  in  the  limit,  of  convergence,  but  is  also  maintained  at  every  time 
step.  Because  of  the  large  linear  system  involved,  extensive  computina  capacity  is 
required.  Example  applications  of  the  method  to  the  driven  cavity  problem  and  to 
the  flow  about  a  circular  cylinder  have  indicated  the  utility  of  the  method,  and 
results  of  these  calculations  have  compared  well  with  previous  numerical  solutions 
and  with  experiment.. 


iii 


AFWAL-TR-87-3051 


FOREWORD 

This  report  Is  the  result  of  research  performed  in  the  Computational 
Aerodynamics  firoup.  Aerodynamic  and  Airframe  Branch,  Aeromechanics  Division,  Flight 
Dynamics  Laboratory,  Air  Force  Wright  Aeronautical  Laboratories,  Wright-Patterson 
Air  Force  Base,  Ohio  45433-6553  hv  Dr  Donald  P.  Rizzetta  from  November  1985  to 
May  1987. 

The  author  is  grateful  to  Dr  M.  R.  Visbal  for  several  helpful  conversations. 


i  v 


AFWAL-TR-87-3051 


TABLE  OF  CONTENTS 


SECTION  PAGE 

I  INTRODUCTION  1 

II  GOVERNING  EQUATIONS  2 

III  NUMERICAL  METHOD  5 

IV  STABILITY  ANALYSIS  10 

V  NUMERICAL  EXAMPLES  14 

VI  DISCUSSION  44 

REFERENCES  46 


AFWAL-TR-87-3051 


LIST  OF  FIGURES 

FIGURE  PAGE 

1  Driven  Cavity  Geometry  and  Boundary  Conditions  15 

9  Driven  Cavity  Convergence  History  for  R  =  100  17 

3  Driven  Cavity  Streamfunctlon  Contours  for  R  =  100  18 

4  Driven  Cavity  Vorticlty  Contours  for  R  =  100  19 

5  Driven  Cavity  Convergence  History  for  R  =  1,000  21 

6  Driven  Cavity  Streamfunctlon  Contours  for  R  =  1,000  22 

7  Driven  Cavity  Vorticlty  Contours  for  R  =  1,000  23 

8  Driven  Cavity  Convergence  History  for  R  =  10,000  24 

9  Driven  Cavity  Streamfunctlon  Contours  tor  R  =  10,000  25 

10  Driven  Cavity  Vorticlty  Contours  tor  R  =  10,000  26 

11  Comparison  of  Driven  Cavity  Velocity  Profiles  with  Solution  of 

Ghia,  Ghia,  and  Shin  28 

12  Circular  Cylinder  Computation  Grid  30 

13  Circular  Cylinder  Grid  Structure  at  Leading  Edge  of 

Computational  Domain  31 

14  Circular  Cylinder  Boundary  Conditions  33 

15  Time  History  of  Lift  Coefficient  for  Circular  Cylinder  36 

16  Time  History  of  Drag  Coefficient  for  Circular  Cylinder  37 

17  Instantaneous  St.reamfunction  Contours  for  Circular  Cylinder  40 

18  Instantaneous  Vorticity  Contours  for  Circular  Cylinder  41 

19  Comparison  of  Drag  Coefficient  for  Circular  Cylinder  with 

Experiment  42 

20  Comparison  of  Strouhal  Number  for  Circular  Cylinder  with 

Experiment  43 


AFWAL -TR-87 - 305 1 


SECTION  I 

INTROOIICTTON  \ 

The  two-dimensional  laminar  incompressible  Navier-Stokes  equations  provide  the 
proper  mathematical  description  for  a  wide  class  of  complex  and  practical  fluid  flow 
problems.  Thus,  this  system  has  been  solved  numerically  by  a  variety  of  techniques  j 

to  obtain  solutions  for  both  steady  and  unsteady  flows  F3-8,  10,  12-15,  19]. 

Historically,  the  technical  sophistication  associated  with  these  techniques  has 
evolved  in  correspondence  with  the  capability  of  computing  hardware.  Decreasing 
computer  storage  limitations  and  increasing  central  processing  rates  have  directly  j 

promoted  the  emergence  of  new  algorithms.  With  the  advent  of  large  scale 
vector-processing  machines,  we  find  it  is  now  possible  to  consider  procedures  which 
were  precluded  hy  previous  restrictions. 

I 

The  present  algorithm  consists  of  a  fully  implicit  finite-difference 
approximation  to  the  governing  partial  differential  equations  and  results  in  a  large 
system  of  linear  equations  representing  the  dependent  variables  at  each  grid  point 
in  the  computational  domain.  This  linear  system  is  then  solved  directly  to  either  I 

advance  an  unsteady  solution  in  time  or  to  proceed  to  the  converged  state  of  a 
steady  flow.  Because  of  the  large  number  of  linear  equations  which  must  be  solved 
repeatedly,  this  formulation  previously  was  inappropriate  for  implementation  on 
computing  machines  having  a  limited  capacity.  When  employed  on  the  current  i 

generation  of  large-scale  computers,  however,  it  may  not  be  impractical.  In 
addition,  the  algorithm  possesses  certain  inherent  advantages  over  other  methods. 

For  unsteady  problems,  formal  second-order  temporal  accuracy  may  be  maintained.  The 
method  contains  no  factorization  error  which  is  commonly  incurred  by 
alternating-direction-implicit  methods.  By  employing  a  direct  method  for 
solution  of  the  linear  equations,  the  use  of  relaxation  parameters  is  avoided,  along 
with  the  appearance  of  spurious  time-like  terms  in  the  finite-difference  equations. 

For  the  case  of  steady  flows  the  iteration  procedure  corresponds  directly  to 
Newton's  method  for  solving  systems  of  nonlinear  partial  differential  equations,  and 
therefore  is  capable  of  achieving  rapid  rates  of  convergence. 


1 


AFWAL-TR-87-3051 


section  II 

GOVERNING  EQUATIONS 


For  an  incompressible  fluid  of  constant  properties,  the  nondimensional 
continuity  and  momentum  equations  may  be  written  in  a  two-dimensional  Cartesian 
coordinate  system  as 


3u 

9x 


=  0 


(1) 


9u  ,  1_ 

at  ax 


(u2)  + 


(2) 


+  JL  (uv)  +  JL  (v2)  +  11 

at  3x  '  ‘  ay  '  ’  3y 


a2v 

3y2 


(3) 


Here,  x  and  y  are  the  independent  Cartesian  coordinates  nondimensional ized  by  the 
reference  length  L,  and  u  and  v  are  the  correspondina  velocity  components 
nondimensional ized  by  the  reference  speed  u^.  L/u^  has  nondimensionalized  the  time, 
t,  P  is  one-half  the  local  pressure  coefficient,  and  R  =  uj./y  is  the  Reynolds 
number.  As  is  customary,  defining  a  streamfunction,  if?,  can  satisfy  the  continuity 
equation  whereby 

«■£  <«> 


v 


,li 

ax 


(5) 


The  pressure  is  eliminated  by  cross-differentiation  of  Eqs  (?)  and  13),  and  the 
momentum  equations  may  be  expressed  in  terms  of  the  vorticity,  w,  defined  as 


w 

I 


a v 
ax 


3u 

ay 


(6) 
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Thus  Eqs  (l)-(3)  are  replaced  by 


which  are  commonly  referred  to  as  the  streamfunction-vorticity  formulation  of  the 
Navier-Stokes  equations.  Equations  (7)  and  (8)  provide  a  concise  representation  of 
the  governing  equations  so  that  any  fluid  flow  field  may  be  defined  by  the  functions 
and  id.  Elimination  of  the  pressure  has  reduced  the  number  of  dependent  variables 
from  three  to  two,  thus  making  the  equation  set  attractive  for  many  numerical 
computations. 

Since  the  early  pioneering  effort  of  Thom  [19],  relaxation  procedures  were 
employed  to  obtain  steady  flow  solutions  for  many  problems  [3,4,151.  Such 
techniques  were  often  slow  to  converge  and  typically  reached  a  limit  beyond  which 
residual  errors  could  no  longer  be  reduced.  These  problems  were  magnified  as  the 
Reynolds  number  increased.  An  effective  means  of  obtaining  iterated  solutions  to 
the  steady  streamfunction-vorticity  equations  has  more  recently  been  achieved  via 
Newton's  method  [5,  6,  131.  By  solving  the  unsteady  equations  of  motion  it  was 
possible  to  attain  steady-state  results  in  the  long-time  limit.  Current  schemes  of 
this  type  can  employ  the  alternating-direction-impl icit  approach  [14]  and  multigrid 
algorithms  [8].  In  the  case  of  nonsteady  flows,  early  methods  utilized  explicit 
procedures  for  time  integration  [7,  10]  which  were  subject  to  stability  limitations. 
The  use  of  upwind-differencing  techniques  aided  in  overcoming  such  restrictions  at 
the  expense  of  degraded  spatial  accuracy.  Contemporary  efficient  algorithms  employ 
an  alternating-direction-implicit  treatment  of  the  vort.icity  transport  equation 
coupled  with  a  direct  solution  of  the  streamfunction  equation.  This  approach, 
however,  is  only  first  order  accurate  in  time.  One  obstacle  associated  with 
obtaining  solutions  to  the  unsteady  equations  of  motion  is  the  lack  of  a  time 
derivative  appearing  in  the  streamfunction  equation,  thus  precluding  a  concise 
representation  which  provides  second-order  temporal  accuracy.  Mehta  and  Lavan  [121 
resolved  this  difficulty  by  solving  Eqs  (7)  and  (8)  simultaneously;  however,  the 
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method  required  iteration  at  every  time  step.  Various  other  procedures  exist  for 

both  the  steady  and  unsteady  equations  of  motion,  but  these  are  too  numerous  to 
mention  here. 
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SECTION  III 
NUMERICAL  METHOD 

To  implement  the  present  numerical  method,  the  Cartesian  coordinates  x  and  y 
are  transformed  to  the  computational  system  £(x,  y),  n(x,  y)  in  the  domain  0  <  £  £  1, 

0  £  n  £  1.  Here,  E  and  n  are  boundary  conforming  to  the  particular  geometry  involved 
to  consider  arbitrary  configurations.  In  addition,  the  computational  grid  is  uniform 
in  E  and  n  (i.e.,  AE  =  constant,  AE  =constant)  so  that  all  spatial  derivatives  may  be 
approximated  by  simple  second-order  accurate  finite-difference  expressions,  as  is 
customary.  Under  this  transformation,  the  governing  equations  are  written  in  conserva¬ 
tive  form  as 


*-  4 

« 


I 


where  E  ■.  5  .  n  .  n  are  the  metric  coefficients  and  J  =  £  n  -  l  n  is  the  Jacobian 
xyxy  xyyx 

of  the  transformation.  For  convenience  we  assume  that  E  and  q  are  orthogonal,  although 
this  simplification  is  not  fundamental  to  the  method.  We  now  presume  that  \p  and  w  are  known 
at  time  level  tn  or  at  the  n-th  iteration  and  Eqs  (9)  and  (10)  are  used  to  advance  the 
solution  to  then  +  1  level.  Approximating  Eq  (9)  accomplishes  this  at  the  n  +  1/2  time 
level  for  the  unsteady  equation  or  at  the  n  +  1  iteration  if  ~  =  0.  For  either 

o  L 

situation,  Eq  (10)  is  satisfied  at  the  n  +  1  level.  Upon  introducing  the  "delta  form" 

of  the  dependent  variables  [2],  i.e.,  Aijj  =  tjjn+^  -  ij/1.  Aw  =  wn+"'  -  ton,  into  the  governing 

2 

equations  and  neglecting  terms  of  order  A  ,  Fqs  (9)  and  (10)  are  expressed  as 


« 


« 
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where  e  =  1  for  unsteady  Flow  and  e=  0  for  steady  flow.  Equations  (11)  and  (1?)  are 
solven  repeatedly  at  each  time  step  or  iteration  to  obtain  Aip  and  Aw,  the  values  of 
and  a)  are  updated,  and  the  procedure  is  repeated  to  advance  the  solution.  The 
unified  formulation  expressed  by  Eqs  (11)  and  (1?)  is  either  second  order  accurate 
in  time  for  nonsteady  flows  or  corresponds  to  Newton's  method  of  iteration  for 
solving  nonlinear  partial  differential  equations  Till  in  the  stationary  case. 

All  spatial  derivatives  appearing  in  Eqs  ^11)  and  (1?!  are  approximated  by 
their  common  centered  second-order  accurate  finite-difference  representations  in 
terms  of  the  dependent  variables  at  each  grid  point  (i,j)of  the  computational  mesh, 
where  for  example  ip.  .  =  ipfc*  ru).  By  defining 

1  9  J  1  9  J 


(2  -  e)  R 

a  ~  4A?An 

(15) 

8  =  1  -  f 

(16) 

R 

Y  "  2A£An 

(17) 
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the  vorticity  transport  equation  becomes 

“  [(“"-l.  3  '  “1,  1-1 1  1-1  +  HI  -  Vl,  j>  ''*1-1,  HI 

+  1-1  ■  “HI,  j)  ''♦hi,  j-1  4  <“H1,  0  -  Hl>  ''♦hi.  Hi] 

+  [“(*?-l.  l-l  -  ♦?-!,  HI 1  -  8(X1.  1  4  XH1,  J>]*Vl.j 

4  [“'♦?+!,  1-1  -  *H1,  1-1 1  -  »  <Y1,  1  4  V  M>]  i“(,  1-1 

4  [In— 4 « (*H, ,  1 4  2X1 ,  J 4  xt_, .  J  * ,  Hi  4  2xi ,  1 4  xi ,  j-1  >]*,.! 

4  [“<♦?-!,  HI  •  ♦hi.  hi  1  -  6  (Yi,  HI  4  yi,j>]  4“1,  Hi 

4  [“<*H1,  HI  -  ♦hi,  1-1  >  -  6  <XH1,  1  4  Xi,  i'J^HI.  1 

■  Y  [‘*1-1.  HI  '  J  •  '♦hi,  1-1  -  *"-l.  1-1  K,  1-1 

+  ^i+1,  j+l  "  ^i-1 ,  j+l  *  “i,  j+l  '  ^i+1,  j+l  "  ^i+l,  j-1  ^  Wi+1,  j 
+  (xi,  j +  xi-i,  j>  “".I,  j +  (yi,  j +  v  j.i) j.i 

^Xi+1,  j  +  2Xi,  j  +  Xi-1,  j  +  Yi,  j+l  +  2Yi,  j  +  Yi,  j-1^  “i,  j 

t(V.,Hl  +V1.3,“?.  HI  +(XH1,  J+x1,  1>“H1,  j] 


(is: 
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Similarly,  the  streamf unction  equation  is 

(X1.JtXl-l.  j  +  (Vi,  J+Y1.Ml4»l.i-l 

(X1*1.  J+XX1.  J  +  X1-l,i  +  Vi,  j*i  +  2  Yi . j  +  ¥i ,  m»*i,  j 

(Vi,  j*l  +  V  J>  j*l  +  <XH1,  j  +  xi,  j>  *i*i.  j  +(3— )^,  j 

<Xf,  j  +  XM.  j>  <*",  j  •  <-l,  J>  -  (XH,,  j  *  Xi,  J>  <*?♦,,  j  -  *?.  j> 

<V  j  +  yi,  j-11  <*!.  j  •  *i.  M>  -  <vi,  j+1  *  Yi,  j>  <*?.  J+i  -  ♦?.  j) 


».r/- j 


(19) 


Equations  (18)  and  (19)  represent  the  linear  system  which  must  be  solved  for 

all  Aip.  .  Au>.  .at  each  time  step.  If  the  computational  domain  consists  of  (T,J) 

1  »j ,  1  *  J 

grid  points  in  ( £,  n)  respectively,  then  the  system,  includinq  boundary  conditions, 

is  comprised  of  2T.1  equations  for  the  variables  Atp.  Aw.  .  where  1  <  i  <_  I  and 

1  ♦  J  1  *  J 

1  <  j  <  J.  This  system  may  be  expressed  notationally  by  the  matrix  equation 

M  I  =  r  (20) 

where  M  is  a  square  matrix  of  dimension  ( 21 J  x  21.1),  and  the  solution  if- and 
right-hand-side  r  are  vectors  of  length  (2IJ).  Although  the  matrix  M  is  large,  it 
is  quite  sparse.  Thus,  we  do  not  find  it  necessary  to  store  a  vast  number  of 
coefficients  to  solve  the  linear  system.  In  addition,  M  is  rich  in  structure  so 
that  a  direct  solution  to  Eq  (201  may  be  obtained  by  efficient  means.  Various  forms 
for  the  matrix  M  are  possible  depending  upon  how  the  equations  are  ordered.  For  the 
moment,  if  we  consider  Dirichlet  conditions  at  all  computational  boundaries  and  let 


3  =  {[(A^  ]  Aw^  1  ). . .  (A^  ^  .  Au)^  j)...(A^  j  Aw^  j)]  [(A^  ,  Au>^  })... 

. . .  (AiPji  ^  j  j)]  C  ( Aipj  ^  1  Ao»j  ^ .,).,.  (A^  .  Ao)j  ^  j)...(A«PIt  j  Auj  ^  j)]J  T  (21) 
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then  M  has  the  block-tridi agonal  *orm 


M  = 


Dt  E, 


°2  E2 


'1-1  UI-1  CI-1 


'I  -J 


(22) 


Here  the  blocks  C.,  n. ,  and  are  (?«1  x  ?»l)  square  matrices  which  are  also  sparse. 

Each  block  itself  has  a  block-tridiagonal  substructure  comprised  of  (2  x  2)  square 

blocks.  Equivalently,  the  C.,  ,  and  E^  matrices  may  be  considered  heptadiagonal . 

Eornberg  [5,  6]  presents  and  alternative  structure  for  the  matrix  M  by  considering  a 

-► 

different  ordering  of  the  solution  vector  1  . 


The  linear  system  represented  by  Eq  (20)  is  solved  by  a  direct  method.  A 
noniterative  technique  is  essential  and  should  be  free  from  limitations  of 
convergence  as  the  mesh  size  is  refined.  In  addition,  a  high  order  of  accuracy  can 
be  obtained  in  the  solution  at  each  time  step,  and  spurious  time-like  relaxation 
errors  will  be  precluded.  We  emphasize  that  the  specific  technique  employed  to 
solve  the  linear  system  is  not  basic  to  the  numerical  method.  Many  efficient 
solution  procedures  are  available  for  this  purpose.  For  the  present  results,  which 
serve  only  as  illustrative  examples,  a  block-Gaussian  elimination  [91  was 
implemented  where  C.,  D.,  and  E.  were  considered  as  full  matrices  without  further 
simp! ification. 
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SECTION  IV 


STABILITY  ANALYSTS 


For  the  steady  case  (e  =  0),  because  the  algorithm  corresponds  to  Newton's 
method,  convergence  of  the  solution  Is  guaranteed  provided  the  Initial  iterate  is 
sufficiently  close  to  the  desired  root.  Mathematical  proof  of  this  property  exists 
TP"1  so  that  no  further  analysis  of  convergence  is  required.  In  addition,  since 
Newton  iteration  is  a  second-order  method,  the  approach  to  the  steady  state 
is  quadratic,  i.e.,  A^n+*  =  0  (Aipn)^.  For  the  unsteady  case  (e  =  1 1 ,  stability  of 
the  procedure  may  be  investigated  by  considering  a  linear  model  system  of  equations 
with  constant  coefficients.  While  the  stability  of  linear  constant,  coefficient, 
systems  is  neither  a  necessary  nor  sufficient,  condition  for  the  stability  of 
corresponding  nonlinear  systems,  it  can  serve  as  a  guide  for  stability  in  the 
nonlinear  situation. 


Equations  (7)  and  (8)  are  modeled  by  the  linear  system 


9w  .  9^  ,  a  9^  ,  l  9ui  .  l  9u>  _  1  ( 9  t 

at  1  ax  2  a7  *  bi  n  *  b2  ay  "  id  ^  + 


a2 
9  0) 

9y2 


dx  dy 


(23) 


(24) 


where  aj,  a,,,  bj,  b?,  are  real  constant  coefficients.  The  finite-difference 
representations  of  these  equations,  which  corresponds  to  the  method  described  in  the 
previous  section,  are  expressed  notationally  as 


n+l  -  n 

*  <a,  5x  +  a2 


(i  n+l  ,  ,n  \  /  n+l  .  n 

5x^2  V(^“ 


.  ,  n+l  ,  n 

R  (6xx  +  6yy)(  ~  2 


(6  +  6  )  +  o)n+1  =  0 

v  xx  yy 


(25) 

(26) 
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Here  6  ,  5  ,  6  ,  6  ,  are  the  common  centered  second-order  accurate  finite- 

a  y  yy 

difference  operators.  Because  the  system  is  linear,  the  "delta"  form  of  the 
equations  will  be  identical  to  Eqs  (25)  and  (2fi)  but  has  been  omitted  for 
convenience.  Upon  defining 

a-.  At  a2  At  b.  At  b?  At  _  At 

al  =  ~ 2~  •  a2  =  ~T~  ’  B1  =  ~2~’  B2  b  2  ’  Y  2R  (27) 

the  system  may  be  written  as 

((*1  Sx  +  a2  6y)  <j>n+1  +  [1  +  (B1  6X  +  B2  6y)  -  Y  (6xx  +  6yy)]  u>n+1 
=  -  (c^  6X  +  a2  6y)  /  +  [1  -  (8]  6X  +  B2  6y)  +  Y  (6xx  +  6yy)]  to"  (28) 

{6xx  +  6yy)  ^+1  +  “n+1  =  0  (29) 


where  we  note  that  Y  :>  0. 

The  algorithm  represented  by  Eqs  (28)  and  (2Q)  is  now  analyzed  for  stability  by 
the  standard  von  Neuman  technique.  We  assume  that  the  difference  equations  are  to 
be  solved  on  a  finite  computational  domain  with  the  uniform  grid  spacings  ax.  Ay. 
Owing  to  the  fact  that  the  system  is  linear,  ip  and  to  may  be  expressed  as  complex 
Fourier  series  of  the  form 


V  y,  t) 


z 

k=0 


Z 

1=0 


i(ekx  +  y ) 


(30) 


w(x,  y,  t)  =  Z  Z  B.,(t)  i(V  +  hy) 
k=0  1=0  Kl  e 


(31) 
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where  i  =  /T,  Afcl  and  Bk1 

for  the  region  0  £  x  £  , 

to  the  following  result: 


are  the  Fourier  coefficients,  and  9k 
0  £  y  £  !_£.  Substitution  of  (30)  and 


*1 


(31)  into  (2ft)  leads 


% 


i Ax  is  in  ^0k  +  VAv/Jcin  (4> 


Ay /sin 


>i  4),)]i} 


4-1  Ay) 


1  i(ekx  +  4>iy) 
kie 


)Aj,e  «(efcx  +  y) 


Ay  Is  i  n 


'<f>i  Ay^ 
,  2 


'f] 

LAxisin 


in  (0k  4x1  +  (s^sin  <♦,  Ay)  ,  j  B"kle  '(9kX  +  *,y> 


(32) 


Due  to  the  orthogonality  of  the  eigenfunctions  in  the  series,  Eq  (32)  can  be 
satisfied  only  if  the  coefficient  of  each  eigenfunction  vanishes  identically, 
so  that  the  coefficients  are  related  according  to  the  expression 


V  Akl  +  ^  +  r°  +  Bki  =  “  VAkl  +  ^  '  Y°  -  V  )  Bkl 

where  ,  u2»  and  a  are  real  and  a>  0,  Similarly,  if  Eqs  (30)  and  (31)  are 
substituted  into  Eq  (29),  and  orthogonality  is  invoked,  there  results 


or 


T  2  2 
-  (— )  2 

/Ax'si  n 


sin 


+  B 


n+1 

kl 


=  0 


a 


(34) 

(35) 


12 


AFWAL-TR-87-3051 


which  also  implies 


Substitution  of  Eqs  (35)  and  (36)  into  Eo  (33)  leads  to 

Ml i  A^1  +  (1  +  yo  +  M2i)a  A^1  =  -vy  +  (1  -  y a  -  jy  )a  A^ 

Equation  (37)  is  of  the  form 

[a(l  +  ya)  +  Ai  ]  AAJT  =  [a(l  -  yo)  -  Ai  ]  A^ 

where  A  is  real.  The  amplification  factor,  g,  is  then  given  by 


(36) 


(37) 


(38) 


An+1 

Akl 

O'  (1  -  y o)^  + 

An 

Akl 

9  $ 

a  (1  +  yo)  + 

1/2 


(39) 


Recalling  that  y,  o  >  0,it  follows  that,  g  <  1  for  all  Ax,  Ay,  At,  and  R  which  is  the 
criterion  for  stability.  Note  that  (35),  (36),  and  (39)  also  imply 


,n+l 

3kl 


3k1 


(40) 
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SECTION  V 


NUMERICAL  EXAMPLES 


To  demonstrate  the  application  of  the  present  algorithm  to  the  solution  of 
physical  problems,  several  numerical  examples  are  computed.  These  include  the 
classic  driven  cavity  problem  and  the  flow  about  a  circular  cylinder,  which  exhibits 
unsteady  vortex  shedding.  Selection  of  these  examples  was  motivated  by  a  desire  to 
illustrate  the  computation  of  both  steady  and  unsteady  flow  fields  and  to  provide 
comparison  with  previous  calculations  and  with  experiment.  Note  that,  for  all  of  the 
results  presented  here  the  metric  coefficients  were  evaluated  numerically  by 
second-order  accurate  formulae  which  correspond  to  centered  differences  at.  all 
interior  grid  nodes  and  three-point  one-sided  differences  at  computational 
boundaries.  By  employing  the  relationships 


J  = 


xSyo 


(41) 


?x  =  yn  J’  =  _x"  J’  nx  =  'V*  ny  =  X£J  (42) 

these  differences  were  carried  out  on  the  computational  grid  (£,  n)  as  is 
conventional . 


Driven  Cavity 


Since  the  first  solutions  were  produced  by  Burggraf  T3l,  flow  in  a  driven 
cavity  has  served  as  an  example  for  the  validation  of  many  algorithms.  The  reasons 
that  this  problem  is  attractive  are  that  the  geometry  is  simple,  a  Cartesian  grid 
system  is  optimal,  and  the  resulting  solutions  exemplify  nonlinear  fluid  behavior, 
including  reversed  flow  regions.  Figure  1  depicts  a  square  cavity  with  sides  of 
length  L  and  with  t.he  origin  located  at  the  lower  left-hand  corner.  The  upper 
boundary  is  moving  from  left  to  right  wifh  the  reference  speed  Boundary 
conditions  on  the  sol  id  walls  consist  of  the  definition  of  the  reference  streamlinp 
and  the  no  slip  condition  for  the  velocity  components.  The  corresponding  boundary 
conditions  for  vorticity  are  in  fact  the  compatabil ity  relations  obtained  bv 
employing  the  limiting  values  of  the  velocity  components  in  the  vorticity  definition 
(Eq  (12))  at  the  walls.  A  second-order  accurate  spatial  representation  of  these 
conditions  results  in  four-point  one-sided  differences.  While  implementation  of  the 
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wall  expressions  poses  no  particular  problem  in  the  n-direction  ( since  the  [h, 
and  E.j  matrices  are  considered  dense),  it  does  alter  the  structure  in  the 
^-direction  and  results  in  the  modified  form  for  the  matrix  M, 


D1  E1  F1  G1 


C2  °2  E2 


M  = 


CI-1  °I-1  EI-1 


AI  BI  CI  DI  — I 


(43) 


Equation  (20)  may  still  be  solved  by  a  block-Gaussian  elimination,  which  for  the 
form  of  M  given  by  Eq  (43)  is  only  slightly  more  complex  than  that  given  by  the 
basic  form  of  Eq  (221.  V/e  consider  it  essential  to  the  method  to  treat  the  wall 
boundary  condition  implicitly.  Only  in  this  way  is  second-order  temporal  accuracy 
maintained  for  unsteady  problems,  or  Newton's  method  implemented  exactly  in  the 
steady  case. 


All  of  the  results  presented  here  pertaining  to  the  driven  cavity  problem  were 
generated  on  uniform  (51  x  51)  computational  grid  with  ax  =  Ay  =  A£  =  An  =  0.02. 
Figure  2  displays  the  convergence  history  for  a  Reynolds  number  of  100  employing 
Newton  iteration  to  attain  the  steady  state  where  initial  conditions  were  taken  as 
homogeneous.  The  residual  is  represented  in  terms  of  the  L ?  norm,  which  is 
evaluated  on  the  unit  computational  square  as 


norm 


I  J 
L  Z 
i=l  j=l 


(44) 


He  note  that  the  limit  of  machine  accuracy  is  achieved  in  approximately  eight 
iterations,  although  plotting  accuracy  may  be  obtained  in  about  four  steps.  Figure  3 
presents  streamfunction  contours  for  this  case  indicating  the  vortex  core  near  x  = 
0.60,  y  =  0.75  as  well  as  the  regions  of  secondary  flow  in  the  lower  corners. 
Corresponding  contours  of  vort.icity  are  shown  in  Fig.  4  and  compare  favorably  with 
other  numerical  computations  f3,  8]. 
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As  indicated  in  Fig.  2,  the  Newton  method  is  always  capable  of  attaining 
a  quadratic  rate  of  convergence  if  the  initial  solution  is  sufficiently  close  to  the 
final  state.  For  arbitrary  fluid  flow  problems  it  is  not  usually  possible  to 
determine  initial  solutions  which  satisfy  this  condition  a  priori.  Tn  the  case  nf 
the  driven  cavity  with  R  =  1,000,  the  Newton  method  did  not  converge  when 
homogeneous  initial  conditions  were  employed.  If,  however,  the  solution  was  allowed 
to  evolve  from  the  initial  homogeneous  state  by  time  integration  for  50  time  steps 
with  At  =  0.5,  machine  accuracy  could  then  be  attained  in  an  additional  seven  steps 
of  Newton  iteration.  The  procedure  of  time  integration  augmented  by  iteration 
illustrates  the  utility  of  the  unified  algorithm.  As  an  alternative  approach, 
Fornberg  T5,  61  computed  the  steady  symmetric  flow  about  a  circular  cylinder 
obtained  via  Newton's  method  by  generating  a  series  of  solutions  with  increasing 
Reynolds  number,  where  each  new  converged  result  served  as  the  initial  solution  for 
the  next  higher  value  of  R. 

The  convergence  history  for  R  =  1,000  appears  in  Fig.  5.  Streamfunction 
contours  for  this  case  are  displayed  in  Fig.  6.  Note  that  the  vortex  core  is 
located  more  toward  the  center  of  the  cavity  than  was  the  case  for  R  =  100.  Tn 
addition  the  regions  of  secondary  flow  in  the  lower  corners  are  more  extensive. 
Although  not  discernable  in  the  figure,  a  secondary  flow  region  is  also  present  in 
the  upper  left-hand  corner.  Corresponding  contours  of  the  vorticity  are  shown  in 
Fig.  7. 

A  more  severe  test  of  the  method  presented  here  is  the  driven  cavity  problem 
for  a  Reynolds  number  of  10,0(10.  In  this  case,  the  converged  result  for  R  =  1,000 
was  employed  as  the  initial  profile  and  the  solution  was  evolved  in  the  time 
accurate  mode  for  250  steps  with  At  =  0.1.  At  this  point  the  solution  was 
sufficiently  close  to  the  final  state  so  that  convergence  to  machine  accuracy  was 
obtained  with  less  than  an  additional  15  steps  of  Newton  iteration.  The  convergence 
history  of  this  calculation  is  presented  in  Fig.  8.  Resulting  streamfunction 
contours  are  shown  in  Fig.  9  and  indicate  that  the  vortex  core  has  continued  to  move 
toward  the  center  of  the  cavity  as  the  Reynolds  number  is  increased.  He  also  note 
that  the  lower  right-hand  secondary  region  has  decreased  in  extent,  while  those  in 
the  upper  and  lower  left-hand  corners  have  increased.  The  corresponding  vorticity 
contours  shown  in  Fig.  If)  indicate  regions  of  extremely  high-vorticity  gradient 
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along  the  upper  and  right-hand  walls  while  the  center  of  the  cavity  is  dominated  by 
an  essentially  constant  vorticity.  While  the  physical  behavior  described  here  is  in 
qualitative  agreement  with  previous  calculations  [8],  the  (51  x  51)  computational 
grid  used  for  this  high  Reynolds  number  flow  is  inadequate  to  resolve  details  with 
much  acuity. 

Comparisons  of  the  present  results  for  all  three  Reynolds  numbers  with  the 
accurate  solutions  of  Ghia,  Ghia,  and  Shin  [8]  are  displayed  in  Fig.  11  in  terms  of 
velocity  profiles  through  the  center  of  the  cavity.  The  solutions  of  Ghia  et  al 
employed  uniform  grids  of  (129  x  129)  for  R  =  100  and  R  =  1,000,  and  a  grid  of  (257 
x  257)  for  the  R  =  10,000  solution.  As  the  Reynolds  number  increases,  the  formation 
of  strong  shear  layers  along  the  upper  and  lower  walls  is  apparent.  Agreement 
between  the  results  is  reasonable  for  the  lower  Reynolds  number  solutions,  but  for  R 
10,000  a  lack  of  resolution  in  the  wall  boundary  layers  of  the  present  calculation 
has  resulted  in  the  indicated  disparity. 

Although  no  numerical  anomalies  were  encountered  for  the  high  Reynolds  number 
case  presented  here  (R  *  10,0001,  this  result  must  be  regarded  as  preliminary  due  to 
the  course  grid  which  was  employed.  Many  algorithms  meet  with  difficulties  as  the 
mesh  size  is  refined.  Based  upon  the  stability  of  the  time  accurate  technique  for 
linear  equations  with  constant  coefficients,  and  the  convergence  properties  of 
Newton's  method,  we  believe  that  spatially  accurate  solutions  can  be  obtained.  Such 
calculations,  however,  will  be  resource  intensive  due  to  the  large  linear  system 
involved. 

Circular  Cylinder 

The  flow  about  a  circular  cylinder  is  another  classic  example  that  has  been 
computed  by  a  number  of  researchers  [4-6,  10,  13],  who  have  obtained  both  steady  and 
unsteady  solutions.  Physically,  steady  flow  fields  are  not  commonly  observed  for 
Reynolds  numbers  above  a  value  of  approximately  40.  Numerically,  however,  it  is 
possible  to  obtain  steady  symmetric  flow  for  much  higher  values  of  R.  The  reason 
for  th-'s  is  that  numeral  freestream  boundary  conditions  and  surface  geometries  may 
be  smooth  to  machine  accuracy.  If  the  differencing  scheme  is  perfectly  symmetric 
fas  in  the  present  easel,  then  there  is  no  possibility  of  asymmetry  evolving  in  the 
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solution.  Unsteady  vortex  shedding  may  be  induced,  however,  by  introducing  a  small 
perturbation  which  will  eventually  form  a  self-sustained  periodic  motion  for  R  >  40. 
This  technique  has  been  employed  in  previous  calculations  [10]  and  was  also 
incorporated  in  the  results  presented  here.  In  the  physical  situation,  there  is 
never  perfect  symmetry  in  either  freestream  conditions  or  the  surface  geometry  so 
that  the  inherently  unstable  flow  will  naturally  exhibit  large  scale  unsteady 
motion. 

The  geometry  to  be  considered  consists  of  a  circular  cylinder  of  diameter  L 
immersed  in  a  freestream  moving  left  to  right  with  the  reference  speed  at  a 
Reynolds  number  of  100.  The  computational  0-grid  is  depicted  in  Fig.  12  and 
consists  of  (121  x  100)  points  in  (£,  n)  where  the  £-di rection  is  circumferential 
around  the  body,  the  n-direction  is  normal  to  it,  and  the  origin  is  located  at  the 
upstream  leading  edge.  A  value  of  Ar  =  0.02618  is  employed  at  the  surface  and  the 
grid  is  stretched  exponentially  to  an  outer  boundary  r  =  100,  r  being  the  radius. 
Circumferentially,  a  uniform  spacing  of  A0  =  3°  was  employed  so  that  the  arc  length 
along  the  surface  was  approximately  equal  to  the  radial  spacing,  i.e.,  r&e  c  Ar. 
These  grid  spacings,  boundary  locations,  and  number  of  points  are  similar  to  those 
employed  by  Jordan  and  Fromm  [10].  The  only  exception  to  a  constant  A0  occurred  at 
the  upstream  leading  edge  whpre  the  first  and  last  mesh  increments  were  adjusted 
so  that  the  S-grid  line  i  *  I  was  located  between  the  lines  i  =  1  and  i  =  2  and  the 
line  i  =  1  was  located  between  the  lines  i  =  1-1  and  1*1. 

The  overlapping  grid  structure  at  the  upstream  outer  boundary  is  illustrated  in 
Fig.  13  for  clarity.  This  construction  results  in  ease  of  application  of  periodic 
boundary  conditions  in  the  circumferential  (£)  direction.  Second-order  accurate 
expressions  are  obtained  by  the  simple  relationships 


♦i.  s  ' 

*1,J  =  (*I.  j* 


l1- 


j  <  J 


(45) 
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with  similar  representations  f0r  w.  These  periodic  conditions  have  once  aqain 
altered  the  basic  structure  of  the  matrix  M  which  now  has  the  form 


M  = 


Periodic  matrices  of  the  type  given  by  Eq  (46)  are  common  in  numerical  problems. 
Several  methods  exist  for  solution  of  the  linear  system  Eq  (?0)  when  M  has  the  form 
given  by  Eq  (46)  ("1,  17,  18"!.  Unfortunately,  all  such  methods  require  considerable 
computational  effort.  Me  estimated  that  computing  time  for  solution  of  the  linear 
system  using  the  periodic  form  of  M  represented  by  Eq  (46)  would  be  approximately 
twice  that  of  the  block-tridiagonal  system  with  M  given  by  Eq  (??).  Furthermore  the 
upstream  region  of  the  flow  field  should  not  experience  any  large  amplitude 
unsteadiness.  For  these  reasons,  the  periodic  conditions  have  been  applied 
explicitly,  which  of  course  are  only  first  order  accurate  in  time.  Values  of  b  and  w 
were  held  constant  along  the  ^-boundaries  over  a  time  step  or  iteration.  After  all 
interior  points  had  been  calculated  at  the  new  level,  boundary  values  were  updated 
according  to 


*1,4*  + 

(2*2,3t*l-M)/3 


1  <  j  <  J 


with  similar  expressions  for  w.  Equation  (47)  follows  directly  from  (46). 


Boundary  conditions  in  the  n-direction  are  indicated  in  Fig  14.  On  the  solid 
surface,  the  no-slip  conditions  were  invoked  and  the  reference  streamline  defined. 
Along  the  forward,  upper,  and  lower  portions  of  the  outer  boundary,  freestream 
conditions  prevailed.  At  the  rearward  outer  boundary,  defined  by  an  arc  of  72° 
bisected  by  the  centerline,  the  conditions  =  u  =  0  were  employed.  The  angle 

n  X 

defining  the  downstream  boundary  was  chosen  rather  arbitrarily.  Other  downstream 
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conditions  were  tested,  namely  =  0  and  second-order  extrapolation  from  interior 

ori 

mesh  points,  and  produced  no  apparent  differences  in  the  computed  results.  We 


should  mention  that  the  condition  =  0  was  applied  with  second-order  accuracy. 

dX 

Admittedly  the  outer  boundary  conditions  are  not  physically  correct.  However,  the 
location  of  this  boundary  is  sufficiently  far  from  the  body  surface  so  that  such 
inaccuracy  did  not  contaminate  the  near  field.  No  numerical  difficulties  resulted 
from  these  conditions,  and  the  appropriateness  of  their  application  can  be  fudged  a 
posteriori  from  the  computed  results. 


For  purposes  of  comparison,  we  find  it  useful  to  generate  the  pressure 
distribution  about  the  body  surface.  This  can  be  obtained  from  the  momentum  Fqs  (2) 
and  (3)  which  were  transformed  to  the  computational  coordinates  (S  ,  nl.  Followed 
the  work  of  Jordan  and  Fromm  [101,  we  see  it  is  convenient  to  formulate  these 
equations  in  terms  of  total  pressure  to  arrive  at 


af  C?  (u?  +  v2)  +  p]  =  1  (nx  -  ny  f£)  +  Jj  (nx  u  +  y 


—  [i  (u2  +  V2)  +  P]  =  -  (f  —  -  F  — )  1 

an  l2  n  j  uy  at  S  at' 


j  (£xu  +  ?yv)»  *  * 


(48) 

(49) 


In  deriving  (48)  and  (49),  note  that  continuity  and  the  definition  of  vorticity  have 
been  employed,  as  well  as  the  orthogonality  of  S  and  n .  At  every  time  step  the 
right-hand  sides  of  (48)  and  (49)  are  available  so  that  the  n-momentum  equation 
may  be  integrated  from  the  outer  boundary,  where  v  =  P  =  0,  ani  u  =  1 ,  inward  along 
y  =  0  to  the  leading  edge.  On  the  solid  boundary  u  =  v  =  0  so  that  (48)  may  be 
int.eorated  from  the  leading  edge  in  the  circumferential  direction  to  obtain  the 
surface  pressure  distribution.  Note  that  along  this  path  Eq  (48)  reduces  to 


3P  1  (ilA)  9s> 

as  "  "  R  \  J  /an 


(50) 
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A  simple  trapezoidal  rule  was  used  to  perform  all  integrations.  The  temporal 
derivatives  appearing  in  Eq  (49)  were  evaluated  using  three-point,  one-sided 
differences.  Having  obtained  the  surface  pressure,  the  lift  and  drag  coefficients 
were  calculated  from  the  following  expressions 


CL  ■  2*  till  <nyP  -  nx  »/R)]|  n  . 

Cn  =  [jl|  (n  P  +  o  t*j/R) ] | 

J  y  n  =  0 


(51) 

(52) 


where  the  indicated  contour  integrations  around  the  solid  body  were  computed  bv 
trapezoidal  rule. 


Homogeneous  initial  conditions  were  prescribed  for  the  dependent  variables  at 
all  interior  mesh  points  and  the  governing  equations  were  integrated  in  time  for  10(1 
steps  with  At  =  0.5.  This  was  followed  by  30  steps  of  Newton  iteration  which  were 
more  than  sufficient  to  obtain  a  steady  symmetric  flow  field  which  was  converged  to 
machine  accuracy.  The  rate  of  convergence  appeared  to  be  slowed  only  slightly  by 
the  use  of  explicit  periodic  boundary  conditions.  Following  the  technique  employed 
by  Jordan  and  Fromm  [10],  unsteady  motion  was  initiated  by  rotating  the  cylinder 
with  a  tangential  velocity  equal  to  20  percent  of  the  freestream  value,  first  in  a 
clockwise  direction  for  10  time  steps  (At  =  0.5)  and  then  counterclockwise  for  10 
time  steps.  At  this  point  the  -forced  rotation  was  terminated  and  integration  in 
time  was  continued  for  an  additional  100  steps  with  an  increment  At  =  0.5.  This 
allowed  the  solution  to  rapidly  evolve  to  a  periodic  self-sustained  unsteady  flow. 
The  time  step  was  then  reduced  to  a  value  of  0.05  to  resolve  the  temporal  content  of 
the  solution.  Other  types  of  initial  perturbation  were  attempted  and  all  appeared 
to  work  equally  well.  The  forced  rotation  employed  here  was  simply  implemented  by 
allowing  a  nonzero  velocity  at  the  cylinder  surface. 

Time  histories  of  the  instantaneous  lift  and  drag  coefficients  for  the 
calculation  appear  in  Fig.  15  and  16  respectively.  Note  that  the  Newton  iterations 
have  been  indicated  by  assigning  a  value  of  At  =  0.5  solely  for  plotting  purposes 
(50  <  t  ^65).  In  the  case  of  steady  symmetric  flow.  Table  1  provides  comparison  of 
the  present  result  with  several  previous  calculations  in  terms  of  the  drag 
coefficient  and  pressure  coefficient  at  the  forward  stagnation  point.  While  the 
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TABLE  1.  Results 

for  Steady  Symmetric 

Flow,  R  =  100 

cn 

Forward  Stagnation  Cp 

Dennis  A  Chang  [4] 

1.056 

1.060 

Jordan  ft  Fromm  [10] 

1.1 

Fornberg  f5,  6] 

1.058,  1.060 

1.061,  1.065 

present  result 

1.044 

1.057 

current  values  are  slightly  lower  than  those  reported  elsewhere,  the  following 
observation  is  significant.  After  the  first  100  time  steps  the  L?  norm  had  dropped 
by  approximately  eight  orders  of  magnitude.  This  amount  of  reduction  would  usually 
be  sufficient  to  satisfy  most  common  convergence  criteria.  At  this  point  the  drag 
coefficient  had  a  value  of  1.058,  which  was  appreciably  lowered  bv  subsequent  Newton 
iteration  (see  Fig.  16,  t  ~  50). 


The  time  histories  of  the  force  coefficients  indicate  that  for  the  unsteady 
flow  which  evolved,  the  solution  was  quite  periodic.  Table  ?  compares  the  present 
calculation  with  the  results  of  Jordan  and  Fromm  HOl  which  employed  a  slightly 
smaller  value  of  the  time  increment,  (At  =  0.015625)  but  had  only  first-order 
temporal  acuity.  Here  the  mean  drag  is  the  average  of  the  maximum  and  minimum 
values,  the  Strouhal  number  is  based  upon  the  period  of  the  lift  cycle,  and  the 
amplitudes  of  the  force  coefficients  are  peak-to-peak  values.  The  time  increment 
employed  in  the  present  computation  (At  =  0.05)  resulted  in  approximately  12 0  time 
steps  describing  one  period  of  the  lift  cycle.  Differences  in  the  amplitudes  of  the 
force  coefficients  between  the  two  solutions,  indicated  in  Table  2,  is  attributed  to 
the  alternative  numerical  methods  employed  in  the  respective  calculations. 
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TABLE  2.  Results  fnr  Unsteady  Flow,  R  =  100 


mean  Cp 

Strouhal  number 

C^  amplitude 

Cp  amplitude 

Jordan  A  Fromm  [10] 

1.26 

0.16 

0.54 

0.012 

present  result 

1.290 

0.162 

0.365 

0.015 

Figure  17  displays  instantaneous  streamfunction  contours  for  the  unsteady  flow 
field.  The  values  of  time  indicated  in  the  figure  are  given  in  radians  and 
referenced  to  an  origin  which  corresponds  to  the  beginning  of  the  last  lift  cycle 
shown  in  Fig.  15.  Evolution  of  the  vortex  shedding  is  apparent.  Antisymmetry 
between  solutions  which  are  tt  radians  apart  in  time  is  also  evident.  Corresponding 
contours  of  vorticit.y  appear  in  Fig.  18.  Although  the  near-body  field  appears  to  be 
adequately  resolved,  the  shed  vorticity  dissipated  in  the  downstream  wake  owing  to 
the  coarseness  of  the  stretched  grid  and  the  homogeneous  downstream  boundary 
conditions.  Comparison  of  the  mean  drag  coefficient  with  the  experimental  data  of 
Wieselsberger  [21]  and  Trit.ton  r?0]  is  provided  in  Fig.  19.  The  present  result  is 
seen  to  fall  well  within  the  experimental  scatter.  In  Fig.  20  the  Strouhal  number 
is  compared  to  the  measurements  of  Roshko  r 161.  We  see  that  the  comparison  is 
adequate  despite  the  fact,  that  the  current  value  may  be  considered  slightly  low. 
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Fiqure  17.  Instantaneous  Streamfunction  Contours  for  Circular  Cylinder 
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Figure  18.  Instantaneous  Vort. 


or 

LU 

o 

QC 

LU 

CD 

CO 

_1 

LU 

CO 

LU 

£ 


O 

E 

or 

<1  o 


3 

CO 

LU 

OC 


LU 

CO 

LU 

q: 

CL 


<b 

- r 

CD 


42 


LOG  |o  (REYNOLDS  NUMBER) 

Figure  19.  Comparison  Drag  Coefficient  for  Circular  Cylinder  with 
Experiment 
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SECTION  VI 

DISCUSSION 

The  numerical  method  presented  here  affords  a  simple  and  straight-forward 
approach  for  obtaining  both  steady  and  unsteady  solutions  to  the  Navier-Stokes 
equations.  In  fact,  it  is  only  the  previous  limitations  of  computing  capacity  which 
have  precluded  its  prior  implementation.  The  unified  formulation  provides  a  means 
of  calculating  time  accurate  solutions  or  of  arriving  at  an  initial  iterate  so 
that  steady  results  may  be  achieved  via  Newton's  method.  Without  this  feature, 
Newton  iteration  may  not  be  useful  for  general  applications  as  there  exists  no 
practical  technique  for  generating  initial  solutions  which  will  converge.  Note  that 
the  method  is  applicable  to  other  sets  of  equations,  especially  those  in 
conservation  law  form,  as  for  example  the  compressible  Navier-Stokes  equations. 

While  the  algorithm  does  possess  certain  attractive  features,  its  use  requires 

extensive  computer  resources.  The  ultimate  efficiency  of  the  method  is  determined 

by  the  solution  technique  employed  for  the  linear  system.  In  this  regard,  the 

block-Gaussian  method  used  here  is  not  considered  to  be  optimal  in  terms  of  either 

computational  time  or  storage.  We  selected  it  because  of  ease  of  impl :~ientat.ion  and 

it  is  not  central  to  the  numerical  algorithm.  Additionally,  it  maintains  certain 

advantages  over  other  solving  procedures.  Because  the  block  matrices  are  considered 

dense,  no  special  logic  is  required  to  satisfy  the  wall  compatibility  relationship 

for  vorticity  in  the  redirect ion.  Although  these  matrices  are  initially  sparse,  the 

storage  they  provide  is  utilized  in  the  solution  process.  All  of  the  results 

presented  here  were  generated  on  a  CRAY  XMP-12  computing  system.  For  the  (51  x  51) 

grid  employed  for  the  driven  cavity  problem  approximately  1. 85-mill  ion  decimal  words 

-3 

of  central  memory  were  required  and  a  data  processing  rate  of  1.0  x  10  central 
processing  unit  (CPU)  second/time  step/grid  point  was  achieved.  Total  computation 
time  to  attain  machine  accuracy  for  each  of  the  cavity  solutions  is  provided  in 
Table  3. 


t 
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TABLE  3. 

Computational  Time  for  Driven 

Cavity 

R  =  ion 

R  =  1,000 

R  =  10,000 

CPU  second 

P.0.8 

148.3 

686.9 

One  final  advantage  of  the  block-elimination  method  is  the  fact  that  external 

storage  devices  easily  may  be  employed  rather  than  central  memory.  Because  only  a 

few  of  the  blocks  need  reside  in  core  at  one  time,  they  may  be  defined,  modified, 

stored  sequentially,  and  later  accessed  sequentially  as  needed.  This  technique  was 

utilized  for  the  cylinder  calculation,  at  the  expense  of  additional  I/O  time.  For 

the  (1?1  x  inn)  grid  employed  for  this  computation,  only  0.58-million  decimal  words 

of  central  memory  were  required  because  of  the  use  of  external  storaqe.  A  data 

3 

processing  rate  of  8.1  x  in  CPU  second/time  step/grid  point  was  attained,  resultinq 
in  a  total  computational  time  of  31.7  CPU  hours  for  the  unsteady  cylinder  flow.  We 
reemphasized  that  the  cited  processing  rates  could  be  considerably  reduced  by  makinq 
use  of  a  more  efficient  solving  technique  than  the  block-elimination  procedure 
employed  here. 
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